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Abstract 


Non-linear entropy stability and a summation-by-parts framework are used to derive entropy stable 
wall boundary conditions for the compressible Navier-Stokes equations. A semi-discrete entropy 
estimate for the entire domain is achieved when the new boundary conditions are coupled with an 
entropy stable discrete interior operator. The data at the boundary are weakly imposed using a 
penalty flux approach and a simultaneous-approximation-term penalty technique. Although dis- 
continuous spectral collocation operators are used herein for the purpose of demonstrating their 
robustness and efficacy, the new boundary conditions are compatible with any diagonal norm 
summation-by-parts spatial operator, including finite element, finite volume, finite difference, dis- 
continuous Galerkin, and ffux reconstruction schemes. The proposed boundary treatment is tested 
for three-dimensional subsonic and supersonic flows. The numerical computations corroborate the 
non-linear stability (entropy stability) and accuracy of the boundary conditions. 
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1 Introduction 


During the last twenty years, scientific computation has become a broadly-used technology in all 
fields of science and engineering due to a million-fold increase in computational power and the 
development of advanced algorithms. However, the great frontier is in the challenge posed by high- 
fidelity simulations of real-world systems, that is, in truly transforming computational science into 
a fully predictive science. Much of scientific computation’s potential remains unexploited-in areas 
such as engineering design, energy assurance, material science, Earth science, medicine, biology, 
security and fundamental science-because the scientific challenges are far too gigantic and complex 
for the current state-of-the-art computational resources [1]. 

In the near future, the transition from petascale to exascale systems will provide an unprece- 
dented opportunity to attack these global challenges using modeling and simulation. However, 
exascale programming models will require a revolutionary approach, rather than the incremental 
approach of previous projects. Rapidly changing high performance computing (HPC) architectures, 
which often include multiple levels of parallelism through heterogeneous architectures, will require 
new paradigms to exploit their full potential. However, the complexity and diversity of issues in 
most of the science community are such that increases in computational power alone will not be 
enough to reach the required goals, and new algorithms, solvers and physical models with better 
mathematical and numerical properties must continue to be developed and integrated into new 
generation supercomputer systems. 

In computational fluid dynamics (CFD), next generation numerical algorithms for use in large 
eddy simulations (LES) and hybrid Reynolds- averaged Navier-Stokes (RANS)-LES simulations will 
undoubtedly rely on efficient high-order accurate formulations (see for instance [2-9]). Although 
high order techniques are well suited for smooth solutions, numerical instabilities may occur if 
the flow contains discontinuities or under-resolved physical features. A variety of mathematically 
stabilization strategies are commonly used to cope with this problem (e.g., filtering [10], weighted 
essentially non-oscillatory (WENO) [11], artificial dissipation, and limiters [2]), but their use for 
practical complex flow applications in realistic geometries is still problematic. 

A very promising and mathematically rigorous alternative is to focus directly on discrete oper- 
ators that are non-linearly stable (entropy stable) for the Navier-Stokes equations. Simultaneously 
satisfying mass, momentum, energy and entropy constraints is a very desirable property. This 
strategy begins by identifying a non-linear neutral stable flux for the Euler equations; a datum 
neither dissipative nor divergent, against which all other operators may be compared. An appro- 
priate amount of dissipation can then be added to achieve the desired smoothness of the solution. 
Regardless of whether dissipation is added, enforcing a senri-discrete entropy constraint enhances 
the stability of the base operator. 

The idea of enforcing entropy stability in numerical operators is quite old [12], and is commonly 
used for low-order operators [13, 14]. An extension of these techniques to include high-order opera- 
tors recently appears in references [15-17] and provides a general procedure for developing entropy 
conservative and entropy stable, diagonal norm summation-by-parts (SBP) operators for the com- 
pressible Navier-Stokes equations. The strong conservation form representation allows them to be 
readily extended to capture shocks via a comparison approach [13, 16]. The generalization to multi- 
domain operators follows immediately using simultaneous- approximation-term (SAT) penalty type 
interface conditions [18], whereas the extension to three-dimensional (3D) curvilinear coordinates is 
obtained by using an appropriate coordinate transformation which satisfies the discrete geometric 
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conservation law [19]. 

Several major hurdles remain, however, on the path towards complete entropy stability of the 
compressible Navier-Stokes equations including shocks. A major obstacle is the need for solid wall 
viscous boundary conditions that preserve the entropy stability property of the interior operator. 
In fact, practical experience indicates that numerical instability frequently originates at solid walls, 
and the interaction of shocks with these physical boundaries is particularly challenging for high 
order formulations. A step towards entropy stable boundary conditions appears in the work of 
Svard and Ozcan [20]. Therein, entropy stable boundary conditions for the Euler equations are 
reported for the far-held and for the Euler no-penetration wall. 

The focus herein is on building non-linearly stable wall boundary conditions for the compress- 
ible Navier-Stokes equations; primarily a task of developing stable wall boundary conditions for 
the viscous terms. At the semi-discrete level, the proposed boundary treatment mimics exactly 
the boundary contribution obtained by applying the entropy stability analysis to the continuous, 
compressible Navier-Stokes equations. Furthermore, the new technique enforces the no-penetration 
Euler wall condition as well as the remaining no-slip and thermal viscous wall conditions in a weak 
sense. The thermal boundary condition is imposed by prescribing the heat entropy how (or heat 
entropy transfer), which is the primary means for exchanging entropy between two thermodynamic 
systems connected by a solid viscous wall. Note that the entropy how at a wall is a quantity 
that in some experiments is directly or indirectly available (e.g., through measurements of the 
wall heat hux and temperature in some supersonic or hypersonic wind tunnel experiments), or can 
be estimate from geometrical parameters and fluid how conditions for the problem at hand. For 
huid-structure interaction simulations (e.g., supersonic and hypersonic how past aerospace vehicles, 
heat-exchangers), the entropy how can be numerically computed at no additional cost while nu- 
merically solving the coupled systems of partial differential equations of the continuum mechanics 
and fluid dynamics. 

Historically, most boundary condition analysis for Navier-Stokes equations is performed at the 
linear level by linearizing about a known state; a rich collection of literature is available [21-23] . The 
non-linear wall boundary conditions developed herein are fundamentally different from those derived 
using linear analysis and cannot rely on a complete mathematical theory. In fact, a fundamental 
shortcoming that limits further development of any entropy stable boundary conditions is the 
incomplete development of the analysis at the continuous level for proving well-posedness of the 
compressible Navier-Stokes equations. Nevertheless, the boundary conditions proposed herein is 
extremely powerful because they provide a mechanism for ensuring the non-linear stability in the 
L 2 norm of the semi-discretized compressible Navier-Stokes equations. In fact, they allow for a 
priori bounds on the entropy function when imposing “solid viscous wall” boundary conditions. 
Furthermore, they are remarkably easy to implement and compatible with any diagonal norm SBP 
spatial operator, including standard finite element (FE), finite volume (FV), finite difference (FD) 
schemes and the more recent class of high-order accurate methods which include the large family 
of discontinuous Galerkin (DG) discretizations [24] and flux reconstruction (FR) schemes [25]. 

The robustness and accuracy of the complete semi-discrete operator (i.e. , the entropy stable in- 
terior operator coupled with new boundary conditions) is demonstrated by computing the subsonic 
and supersonic flow past a 3D square cylinder without the need to introduce artificial dissipation, 
limiting techniques, or filtering to stabilize the computations, a feat unattainable with several al- 
ternative approaches to wall boundary conditions based on linear analysis. In fact, instabilities are 
observed when wall boundary conditions designed with linear analysis are used in combination with 
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highly clustered grids and/or high order polynomials, or extremely coarse grids near solid walls, 
which yield to unresolved physical flow features. 

The paper is organized as follows. In Section 2, the compressible Navier-Stokes equations, their 
entropy function and symmetrized form are introduced. Section 3 presents the entropy analysis of 
the viscous wall boundary conditions at the continuous level. Section 4 provides a discussion of the 
inviscid flux condition which allows the construction of high-order accurate entropy conservative 
and entropy stable fluxes at the semi-discrete level, for the interior operator. The new entropy 
stable wall boundary conditions and their non-linear stability analysis are presented in Section 5. 
Finally, the accuracy and high level robustness of the proposed approach in combination with a 
high-order accurate entropy stable interior operator are demonstrated in Section 6. Conclusions 
are discussed in Section 7. 


2 The compressible Navier-Stokes equations 


Consider a fluid in a domain O with boundary surface denoted by <9P, without radiation and 
external volume forces. In this context, the compressible Navier-Stokes equations, equipped with 
suitable boundary and initial conditions, may be expressed in the form 


a? , vP = 9/T , , 

dt dxi dxi ’ 
q\an = g {B \x,t), xedQ, 
q(x, 0) = g(°\x), x € P, 


, te[o,oo), 
t e [o, oo), 


(i) 


where the Cartesian coordinates, x = (x\,X2-,x 3 )^ S R 3 , and time, t S R, are the independent 
variables. Note that in (1) Einstein notation is used. The vectors q(x,t ) : R 3 x R — > R 4 , f\ 1] = 
and f- V ^ = \q,Vq) are the conserved variables, the inviscid fluxes and viscous fluxes, 
respectively in the i direction. 1 Both boundary conditions, g^ B \ and initial data, g^°\ are assumed 
to be bounded L 2 nL°°. Furthermore, g^ is also assumed to contain (linearly) well-posed Dirichlet 
and or Neumann and or Robin data. The conservative variable vector is 


q = {p,pu 1 ,pu 2 ,pu 3 ,pE) T , (2) 

where p denotes the density, u = (ui,U 2 ,U 3 ) T is the velocity vector, and E is the specific total 
energy, which is the sum of the specific internal energy, e, (defined later) and the kinetic energy, 
\ujUj. The convective fluxes are 


fj; I} = ( pui,puiui + Sap, puiu 2 + S i2 p, puiu 3 + 5 i3 p , pui,H) T , (3) 

where p, H and Sij are the pressure, the specific total enthalpy and the Kronecker delta, respectively. 
The definition of the viscous flux is 

f[ V) = (0, Ta,Ti2, T i3 , TjiUj - q,;) T , (4) 

x The symbol Vg denotes the gradient of the conservative variables. 
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where the shear stress tensor, assuming a zero bulk viscosity, is defined as 


f dui duj 2 duk \ 

\<9ar,- dxi 3 dxk ) ’ 


and the heat flux, according to the Fourier heat conduction law, is 


q; = 



( 5 ) 


( 6 ) 


The symbols T, p = p ( T ) and k = k (T) which appear in (5) and (6) denote the temperature, the 
dynamic viscosity and the thermal conductivity, respectively. The viscous flux terms in (4) can 
also be expressed as 


AV) 


= a 


dq 


V dx 3 


, dv 

C W 


(7) 


where Cij and d- are five- by- five coefficient matrices, which are a function of the solution variables, 

and v = (p 1 u\,U 2 ,U 3 ,T) T is the vector of the primitive variables. The expression of the matrices 
(T is given in Appendix A. As we will show later, relation (7) is a very convenient form for the 
entropy stability analysis. 

The constitutive relations for a calorically perfect gas are 


e — c v T , 


h = H — -Uj Uj = c p T, 


( 8 ) 


where the symbols c v and c p denote the specific heat capacity at constant volume and constant 
pressure, respectively, and 


P = PRT, 


R = 


Rn 


MW ’ 


(9) 


where R u is the universal gas constant and MW is the molecular weight of the gas. The speed of 
sound for a perfect gas is 

c=^RT, ( 10 ) 

Cp ii 

In the entropy analysis that will follow, the definition of the thermodynamic entropy is the explicit 
form, 

s = — rriog (£-) ~ Rl °s (— ), (n) 

7 1 V Too / V Poa ) 

where and are the reference temperature and density, respectively. 

Note that in practical situations, most simulations are performed in computational space, that 
is, by transforming all elements in physical space to standard elements in computational space via 

a smooth mapping. However, to keep the notation as simple as possible, a uniform Cartesian grid is 

considered in the derivations presented herein. The extension to generalized curvilinear coordinates 
and unstructured grids follows immediately if the transformation from computational to physical 
space preserves the semi-discrete geometric conservation (GCL) [19]. 
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2.1 Entropy function and entropy variables of the compressible Navier— Stokes 
equations 

In the linear hyperbolic framework, L 2 stability is sought as a discrete analogue for a priori energy 
estimates available in the differential formulation (e.g., Richtmyer and Morton [26] and Gustafsson, 
Kreiss and Oliger [21]; Kreiss and Lorenz [27]). In the context of non-linear problems dominated 
by non-linear convection, we seek entropy stability as a discrete analogue for the corresponding 
statement in the differential formulation. Moreover, for systems of conservation laws, stability with 
respect to a mathematical entropy function is considered as an admissibility criterion for selecting 
physically relevant weak solutions. In fact, the entropy condition plays a decisive role in the theory 
and numerics of such problems as shown for instance by Lax [28] and Smoller [29]. 

Harten [30] and Tadmor [31] showed that systems of conservation laws are symmetrizable if and 
only if they are equipped with a convex mathematical entropy function. Given a set of conservation 
variables q(x,t), the entropy variables which symmetrize the system are defined as the derivatives 
of the mathematical entropy function with respect to q{x,t). Hughes and co-authors [12] extended 
these ideas to the compressible Navier-Stokes equations (1). There, it is shown that the mathe- 
matical entropy must be an affine function of the physical (or thermodynamic) entropy function 
and that semi-discrete solutions obtained from a weighted residual formulation based on entropy 
variables will respect the second law of thermodynamics. Hence, it is again found that the entropy 
variables are a critical ingredient in the design of numerical schemes exhibiting non-linear stability. 

Definition 2.1. A scalar function S = S(q ) is an entropy function of Equation (1) if it satisfies 
the following conditions: 


• The function S(q) is convex and when differentiated, simultaneously contracts all the inviscid 
spatial fluxes as follows 

ds djf = dsdjfl dq^ = m dq^ = m (12) 

dq dxi dq dq dxi dq dxi dxi ’ ’ ’ 

The components of the contracting vector, dS/dq, are the entropy variables denoted as w T = 
dS/dq. Fflq ), i = 1,2,3 are the entropy fluxes in the three Cartesian directions. 


The new entropy variables, w, symmetrize Equation (1): 


dq , df j I] dfl V) 


dt dxi 


dq dw df- 1 ^ dw d 


dxi 


— "o wr H — 7T — w w — ( Cij — — I — 0, i — 1,2,3 


dw dt dw dxi dxi 


dw \ 


dxj 


(13) 


with the symmetry conditions: dq/dw = (dq/dw) T . df^ /dw = (df^ /dw) andcij = cjj, 

with the matrices Cij included in Appendix A. Because the entropy is convex, the Hessian d 2 S/dq 2 = 
dw/dq is symmetric positive definite (SPD), 


C J 


d 2 S 

dq 2 


C > 0, vc + o, 


(14) 


and yields a one-to-one mapping from conservation variables, q, to entropy variables, w T = dS/dq. 
Likewise, dw/dq is SPD because dq/dw = dw/dq 1 and SPD matrices are invertible. The entropy 
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and corresponding entropy flux are often denoted as entropy-entropy flux pair, ( S,F ), [13]. If the 
entropy function S(q ) is convex, a bound on its global estimate can be converted into an a priori 
estimate on the solution of (1) in suitable L p space [32]. 

The symmetry of the matrices dq/dw , df- 1 /dw indicates that the conservative variables, q, 
and the inviscid fluxes, f^\ are Jacobians of scalar functions with respect to the entropy variables, 


Q 


T 


d(p 

dw' 



dfk 
dw ’ 


(15) 


where the non-linear function, tp, is called the potential and ifi, i = 1,2,3 are called the poten- 
tial fluxes [13]. Just as the entropy function is convex with respect to the conservative variables 
(. d 2 S/dq 2 is SPD), the potential function is convex with respect to the entropy variables. 

Godunov [33] proved that (see reference [30] for a detailed summary of the proof): 


Theorem 2.1. If Equation (1) can be symmetrized by introducing new variables w, and q is a 
convex function of ip (the so-called potential), then an entropy function S(q) is given by 


ip = w T q — S, 


(16) 


and the entropy fluxes satisfy 

A = w T fP~Fi, (17) 

where ipi, i = 1, 2, 3, are the so-called potential fluxes. The potential and the corresponding potential 
flux are usually denoted as potential-potential flux pair, (p,^), [13]. 

In the specific case of the compressible Navier-Stokes equations, the entropy-entropy flux pair 
is 

S=-ps , Fi = —pui s, i = 1,2,3, (18) 

and the potential-potential flux pair is 

p = p R, tpi = puiR, i = 1,2,3. (19) 


Note that the mathematical entropy has the opposite sign of the thermodynamic entropy. To avoid 
confusion, herein entropy refers to the mathematical entropy unless otherwise noted. The entropy 
variables using the pair in (18) are 

dS\ T fh UiUi u\ U 2 us 1\ T 
~dq) ~ \T ~ S ~ ~2T~’ Y' Y' Y'~t) ’ 

and may be shown to have a one-to-one mapping with the conservative variables provided p,T > 0. 
Expressly: 

C T |^C T> °, VC/0, p, T > 0. 

This restriction on density and temperature weakens the entropy proof, making it less than full 
measure of non-linear stability. Another mechanism must be employed to bound p and T away from 
zero to guarantee stability and positivity; positivity preservation will not be considered herein. 




3 No-slip boundary conditions: Continuous analysis 


The problem of well-posed boundary conditions is an essential question in many area of physics. 
For the two- (2D) and three-dimensional (3D) Navier-Stokes equations, the number of boundary 
conditions implying well-posedness can be obtained using the Laplace transform technique [21]; 
a complicated procedure for system of partial differential equations (PDEs) like the compressible 
Navier-Stokes equations. Nordstrom and Svard [22] proposed an alternative semi-discrete approach 
to arrive at the number and type of boundary conditions for a general time-dependent system of 
PDEs. This analysis was applied to the linearized 3D compressible Navier-Stokes equations for 
different flow regimes and the case of the no-penetration velocity condition. 

In 2008, Svard and Nordstrom [23] showed that the no-slip boundary conditions together with a 
boundary condition on the temperature imply stability for the 3D linearized compressible Navier- 
Stokes equations. Their result, can also be generalized to assess the stability of the non-linear 
problem for smooth solutions as indicated in [21,34] and the references therein. In this section we 
address the non-linear stability (entropy stability) of the viscous wall boundary conditions for the 
(non-linear) compressible Navier-Stokes equations given in (1). 

Contracting the system of equations (1) with the entropy variables and using relations (12) and 
(13) results in the differential form of the entropy equation: 


dw 


T 


asdg as at f> as asa/} v > a , T m 

dq dt dq dxi dt d x; dq dxi dxi V 1 ' ' ' l 


( V ) 


d 


,JAV) 




dx i/ 
dw\ T 


dw 


(21) 


dxi J 


dxa 


Integrating Equation 21 over the domain yields a global conservation statement for the entropy in 
the domain 


/■ 


— IS dx = 
dt 

n 


w T fi V) - Fi 


f f dw\ T ^ dw 

.an J \ dxi ) ^ dxj X 

n 


w T fP-F l ] - DT , (22) 

J Oil 


where DT is the viscous dissipation term, 


DT 


f f dw\ T ^ dw 

J \MJ ° ij d ^ 


dx. 


It is shown in [16] that the 15 x 15 matrix c (i.e. , c)j, 1 <i,j < 3) in the integral term is positive 
semi-definite, which makes the term —DT strictly dissipative. Note that the entropy can only 
increase in the domain based on data that convects and diffuses through the boundaries. The sign 
of the entropy change due to viscous dissipation is always negative. 2 

To simplify, we let the domain of interest, D, be the unit cube 0 < x\,X 2 ,xs < 1. Expanding 

2 From a physical point of view, this means that the viscous dissipation always yields an increase of the thermo- 
dynamic entropy. 
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the Einstein notation in equation (22) yields 
d_ 
dt 

+ 


d f 

— S dx i dx 2 dx 3 = — DT 

n 


£ 1=0 


,.T 


+ 


a?i = l 


+ 


/ 


£ 2=0 


+ 


£ 2=1 


+ 


/ 


£3=0 


+ 


£3 = 1 


,,T 


,,T 


,,T 


,,T 


,,T 



dw 



dw 



dw 

Cll 

dx 1 

+ 

Cl 2 

dx 2 

+ 

Cl 3 

dx 3 


dw 

+ 


dw 

+ 


dw 

Cll 

dx 1 

C12 

dx 2 

Cl 3 

dx 3 


dw 



dw 

+ 


dw 

C12 

dx 1 

+ 

C22 

dx 2 

C23 

dx 3 


dw 



dw 

+ 


dw 

C12 

dx 1 

+ 

C22 

dx 2 

C23 

dx 3 


dw 



dw 

+ 


dw 

Cl 3 

dx 1 

+ 

C23 

dx 2 

C33 

dx 3 


dw 



dw 

+ 


dw 

Cl 3 

dx 1 

+ 

C23 

dx 2 

C33 

dx 3 


dx 2 dX3 

cfaq dx3 
cfaq dx3 
dx i dx2 
cfaq dx2 


(23) 


Consider the case of a wall placed at aq = 0, and assume that all the other boundaries terms 
are stable; their contributions are neglected without loosing generality. Therefore, estimate (23) 
reduces to 


d_ 

dt 


j S dx i dx 2 dx 3 = —DT 


n 

Fi 

£ 1=0 



w 


T 


Cll 


dw 

dx\ 


Cl2 


dw 

dx 2 


+ C13 



dx 2 dx 3 . 


(24) 


To bound the time derivative of the entropy, the right-hand-side (RHS) of Equation (24) requires 
boundary data. In the case of a solid viscous wall (assuming linear analysis) four independent 
boundary conditions must be imposed [34] . 3 The no-slip boundary conditions are u\ = U 2 = 113 = 0 
(i.e., the velocity vector relative to the wall is the zero vector). For the linearized compressible 
Navier-Stokes equations, the fourth condition is the gradient of the temperature normal to the 
wall, ( dT/dn) wa ii , (Neumann boundary condition; e.g., the adiabatic wall), or the temperature of 
the wall, T wa u, (the Dirichlet or isothermal wall boundary condition), or a mixture of Dirichlet 
and Neumann conditions (the Robin boundary condition). These four boundary conditions lead 
to energy stability (linear stability); see, for instance, [23,35]. In the remainder of this section, we 
will show the type and the form of the wall boundary conditions that have to be imposed to bound 
estimate (24) and, hence, to attain entropy stability. 

Consider the inviscid contribution to the boundary terms in (24). 


Theorem 3 . 1 . The no-slip boundary conditions u± = U2 = U3 = 0 bound the inviscid contribution 
to the time derivative of the entropy in equation (24) . 

3 A solid wall behaves like a subsonic outflow [22]. 
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Proof. Equation (17) provides the definition of the entropy flux, F\ : 

F\ = w T / x (/) - ip!, tpi = puiR. (25) 

Substituting the no-slip conditions, u\ = U 2 = = 0, into the definition of the inviscid flux, 

(Equation 3) and the condition u\ = 0 into the definition of ipi, yields the desired result F\ = 0. □ 

A similar proof appears in the work of Svard and Ozcan [20], where a SAT approach on the con- 
servative variables is used. 


Consider now the viscous contribution to the boundary terms in (24). 

Theorem 3.2. The entropy flux boundary condition, 

dT 1 

^=Tnr < 26 > 

where the symbol n defines the normal direction to the wall, bounds the viscous contribution to the 
time derivative of the entropy (24) . 


Proof. Using the definition of matrices Cjj given in Appendix A, the viscous vector-matrix-vector 
contraction given in (24) yields the following term 


— w 


^ dw ^ dw ^ dw \ 

C11 w b Ci2w b C13- — 

OX 1 <JX 2 OX 3 J 


K 


dT l\ 
d^TJ 


(27) 


Therefore, specifying the condition g(t) = where g(t) is a known bounded function (i.e., 

L 2 n L°° ), eliminates the last potential source of instability on the right-hand-side of equation 
(24). □ 


The boundary condition given by Theorem 3.2 at first glance appears ad hoc. Note, however, 
that the scalar value k accounts for the change in entropy at the boundary x\ = 0, due to 

the wall heat flux, q wa ii, also denoted heat entropy transfer or heat entropy flow [36]. In fact, 


OT 1 d r 1 

k - — — = — k — — 10(0) — — 
dx 1 T dx 1 w(5) 


= [log(r)1 


Qwall 

T wa ll ’ 


(28) 


where w(5) represents the fifth component of the entropy variable vector, w, in (20). Equation 
(28) indicates that, in the context of the entropy stability analysis of the compressible Navier- 
Stokes equations, it is admissible and physically (thermodynamically) correct to impose the fourth 
wall boundary condition as given in Theorem 3.2. Such a boundary conditions can be seen as a 
non-linear mixed Dirichlet and Neumann boundary condition. 

To the best of our knowledge, Theorem 3.2 provides new insight for future development of any 
entropy stable boundary conditions for the compressible Navier-Stokes equations. 


Remark 3.1. We strongly remark that the non-linear contraction obtained in (27) is different 
from that obtained from the linearized compressible Navier-Stokes equations [23, 35], The linear 
analysis produces velocity gradient terms in the energy estimate, and temperature gradient terms in 
the normal direction of the form Tff^, with T and being perturbations. 

Remark 3.2. The boundary condition g(f) = 0, which corresponds to an adiabatic wall dT/dx 1 = 
0, bounds the solution, and, as physically expected, does not contribute to the time derivative of the 
entropy (24). 
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4 Entropy stable spectral discontinuous collocation method for 
the semi-discretized system: Interior operator 


In this section, we summarize the main results that allow us to construct a numerical high order 
entropy stable discontinuous spectral collocation interior operator of any order, p. The formalism 
provided here, will then be used in Section 5 to design new entropy stable solid wall boundary 
conditions for the semi- discretized compressible Navier-Stokes equations. 

Using an SBP operator and its equivalent telescoping form (see Appendix B), the semi-discrete 
form of the compressible Navier-Stokes equations (1) becomes 


<9q 

dt 


= V, r 



+ CijV Xj q ) + V x . 1 (g {B) 


+ g (7n) ) 


V-}A Z 




+ V~H g (B) + g (Jn) ) 


q(x, 0) = g^(x), iSll. 

(29) 

The source terms g' and g( In '> contain the enforcement of boundary and interface conditions, 
respectively; and g*°) represents the initial condition. The matrix V incorporates the local grid 
spacing into the derivative definition, and it may be thought of as a mass matrix in the context 
of Galerkin FE method. While it is not true in general that V is diagonal, herein the focus is 
exclusively on diagonal norm SBP operators, based on fixed element-based polynomials of order 
p. The matrix T> is used to approximate the first derivative; and it is defined as V~ l Q, where the 
nearly skew-symmetric matrix, Q, is an undivided differencing operator where all rows sum to zero 
and the first and last column sum to —1 and 1, respectively. The operator A is the telescopic flux 
matrix and allows to express the semi-discrete system in a telescopic flux form, by evaluating the 
fluxes at the collocated flux points, and f ( ' \ 4 (see Figure 1). 


fo fl /2 


h 

h 


fz 

Uq 

Ul 


U 2 

II k— 

• 


• 

X\ 

X 2 


X3 

Xq X 1 


X 2 



/3 


* 


X3 


h 

U3 

X4 


fh 


* 


U 4 

II 


X 5 
X 4 X$ 


i ^9 _ rs -i6 

1 10 V 7 45 


0 


+ 16 
45 


+9 

10 


+ 1 


Figure 1. The one-dimensional discretization for p = 4 Legendre collocation. Solution points are 
denoted by • and flux points are denoted by x . 


The semi-discrete entropy estimate is achieved by mimicking term by term the continuous 
estimate given in Equation (21). The non-linear analysis begins by contracting the semi-discrete 
equations given in Equation (29), with the entropy variables: w t T’. For clarity of presentation, 
but without loss of generality, the derivation is simplified to one spatial dimension. Tensor product 
algebra allows the results to extend directly to three dimensions. The resulting global equation 

4 In the reminder of this work, all quantities evaluated at the flux points are denoted by an over-bar. 
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that governs the time derivative of the entropy is given by 


+ w T A fW = w T A f(^) + w T (g (B) + g (/n) )> (30) 

where 

w = (w(qi) T ,w(q 2 ) T , . . . ,w(q N ) T ^j , 
is the vector of entropy variables. 


4.1 Time derivative 

The time derivative in (30) is in mimetic form for diagonal norm SBP operators. The entropy 
variables are defined by the expression w T = dS/dq, which when combined with the definition of 
entropy yields the point-wise expression 

T dqi dSidqi dSi 

w 'm = dim = aP v “' 

Now, define the diagonal matrices dS/dq = W = diag[w]. Since V is a diagonal matrix and 
arbitrary diagonal matrices commute, the semi-discrete rate of change of entropy becomes 


<9q 


<9q 


r <9q 


w'V^ = l'WPv = 1 ' VW^ = 1 T V = l'V=^, 

dt dt dt <9q dt dt 


dS dq 


ib 


.as 


where 

1 = ( 1 , 1 ,..., 1 ) T , 

is a vector with N elements. 5 


4.2 Inviscid terms 

The inviscid portion of Equation (30) is entropy conservative if it satisfies 

w T Af (/) = F(q N ) - F( qi ) = F(q N ) - F( gi ) = 1 T AF. (31) 

Note that in (31), the first and last flux points are coincident with the first and last solution 
points, which enables the endpoint fluxes to be consistent (see Figure 1). One plausible solution 
to Equation (31) is a point-wise relation between solution and flux-point data, which telescopes 
across the domain and produces the entropy fluxes at the boundaries. Tadmor [13] developed such 
a solution based on second-order accurate centered operators. Carpenter and co-authors [37], have 
generalized this solution for Legendre spectral collocation operators of any order of accuracy, p. 

In the following paragraphs, we present, without any proof, the main theorems which allow 
to construct inviscid entropy conservative and stable fluxes of any order of accuracy. Interested 
readers should consult [16, 37] and the references therein for details. Note that in this section the 
subscripts i — 1 , i and i+ 1 are used to denote a scalar or vector quantity at the i — 1 , iori + 1 
collocated point, and do not have to be confused with the subscript used for instance in (1). 

5 N = p + 1 for a pth-order scheme. 
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Theorem 4.1. The local conditions 

(w i+ i - Wi) T Ji = fa +1 - fa, i = 1,2, ...,7V- 1 ; fa = fa, = i]n (32) 

when summed, telescope across the domain and satisfy the entropy conservative condition given 
in Equation (31). A flux that satisfies this condition given in equation (32) is denoted f> . The 
potentials fa+i and fa need not be the point-wise fa+i and fa, respectively. 

A possible strategy for constructing high order entropy conservative fluxes is to construct a 
linear combination of two-point entropy conservative fluxes by using the coefficients in the SBP 
matrix Q. This approach follows immediately from the generalized telescoping structural properties 
of diagonal norm SBP operators given in Section B.l. Because it requires only the existence of 
a two-point entropy conservative flux formula and the coefficients of Q, it is valid for any SBP 
operator that satisfies the constraints given in (B3). Thus, it is also valid for Legendre spectral 
collocation operators used herein. 

The following theorem establishes the accuracy of the new fluxes. 

Theorem 4.2. [16] A two-point entropy conservative flux can be extended to high order with 

formal boundary closures by using the form 

N i 

li S) = J2 J2 2Q ^s(Qe,Qk), l<i<N-l t (33) 

k=i-\-l £=1 

when the two-point non-dissipative function from Tadmor [13], 

i 

Is (■ Qk , Qi) = J 9 (w(qk) + € {w{qe) - w(q k ))) d£, g{w{u)) = f(u), (34) 

o 

is used. The coefficient, Qki, corresponds to the k row and l column in Q. 

Thus, Theorem 4.2 ensures that a high order flux constructed from a linear combination of 
two-point entropy conservative fluxes retains the design order of the original discrete operator for 
any diagonal norm SBP matrix Q. 

The following theorem establishes instead that the linear combination of two-point entropy 
conservative fluxes does preserve the property of entropy stability for any arbitrary diagonal norm 
SBP matrix Q. 

Theorem 4.3. [16] A two-point high order entropy conservative flux satisfying Equation (32) with 
formal boundary closures can be constructed using Equation (33), 

N i 

7f )= E E 2 ^/s(®>%), 1<*<JV-1, (35) 

k=i+l i— 1 

where fs(qe,Qk ) is any two-point non-dissipative function that satisfies the entropy conservation 
condition 

{w£- w k ) T J s (qe,qk) = fa-fa- (36) 
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w T ^- 1 Af (S) = T^AF = ^-F (q) + T p+ 1 , 


T/ie /ng/i order entropy conservative flux satisfies an additional local entropy conservation property, 

(37) 

where T p+ \ is the truncation error of the approximation, or equivalently, 

• T ( ff ' ' - Ji-i) = - ^i-i) , l<i<N, (38) 


w,- 


where 


N i 

Fi = EE I + «+) T fs (%> <?*:) - {A + i’k) 

k=i-\-l 1=1 


1 < i < N - 1. 


(39) 


4.2.1 Affordable entropy consistent Euler flux 


The inviscid terms in the discretization of the compressible Navier-Stokes equations (29) are cal- 
culated according to Equations (33) and (34) by using the two-point entropy conservative flux of 
Ismail and Roe [38], 


fs,j(Qi > Qi+ 1) = ( puj,piijUi + Sjip, piijU 2 + 5j2p, puju 3 + 5j 3 p, phjH^J , 


u = 


Uj , Vj + l 


VT t y/Ti + l 


Pi I Pi+1 

„ FF " r 

— » p = — 


1 > 


FF + 


h = R 


log (vfc) 
1 + 1 


(0i + e 2 ) 


o i = 


+T \/T+i 

VTiPi + V^i+lPi+1 


^2 = 


+ V^+i) y^+i/7+i) 

7 + 1 


log ( x /% 1 


7-1 


log 


H = h + -ugU£, p = 


Ti Pi \ | 1 1 ] 

T+i Pi+i ) \^y/Ti \jTi + 1 y 

+ vfc) “ V^l«+l) 

2 (log(v / T)pi) - log(VT, + i Pm )) 


(40) 


The index j denotes the spatial direction. This somewhat complicated explicit form is the first 
entropy conservative flux for the convective terms with low enough computational cost to be im- 
plemented in a practical simulation code. 


4.2.2 Entropy stable inviscid interface flux 

Herein, the solution between adjoining elements is allowed to be discontinuous. An interface flux 
that preserves the entropy consistency of the interior operators on either side of the interface is 
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needed. An entropy consistent (or entropy conservative) inviscid interface flux constructed accord- 
ing to equations (33) and (34) by using (40) is indicated as f sr (qj \qj + ^. A more dissipative and 

hence stable entropy inviscid interface flux f ssr (q\ \ gj - " 1 J is constructed as 








+ A (™! +) - } ) > 


(41) 


where A is a negative semi-definite interface matrix with zero or negative eigenvalues. The entropy 
stable flux f ssr (q[ \qj + ^ is more dissipative than the entropy conservative inviscid flux, as is 

easily verified by contracting f ssr fqj \q\ + ^j against the entropy variables to yield the expression 




j?ssr 






+ 






(42) 


The superscripts (— ) and (+) combined with the subscript i denote the left and right state used to 
compute the two-point entropy conservative flux and therefore replace the subscripts i and i + 1 in 
(40). The matrix A can be constructed using different approaches, e.g., using an upwind operator 
that dissipates each characteristic wave based on the magnitude of its eigenvalue: 


£»SSC 

d 

dq 

dq 

dw 


,(-) 


Qi 


(+)') _ fsr 


Qi 


(-) „(+) 


+ i/2y\\\y 


T 


(-) (+) 
W „■ — W- 


f (?) = TAT 1 


= TT 


(43) 


where A and T are the diagonal matrix of the eigenvalues and the matrix of the eigenvectors, 
respectively. Note that the relation = TT r is achieved by an appropriate scaling of the 
rotation eigenvectors. See the work of Merriam [39] for more details. Unless otherwise noted, the 
entropy stable characteristic flux is used in all test simulations presented herein. 


4.3 Viscous Terms 

Using the SBP formalism (see Appendix B.2), the contribution of the viscous terms to the semi- 
discrete time derivative of the entropy is 

w T Af^' ^ = w t BciiT ) w — (Pw) T Vc\\ (Dw) . (44) 

The last term is negative semi-definite. As with the continuous estimate given in (23), only the 
boundary term can produce a growth of the entropy, and thus the approximation of the viscous 
terms is entropy stable. (Entropy stable boundary conditions bound these terms.) 


5 Entropy stable solid wall boundary conditions for the semi- 
discrete system 

An estimate for the time derivative of the entropy of a isolated element is derived, followed by a 
derivation of entropy stable penalty terms that imposed physical data on the boundaries. 6 

6 The same boundary conditions (without stability proofs) could be used for almost any spatial discretization 
including the family of DG methods, FR approaches, WENO schemes, FD and FV methods. 
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5.1 General approach for the entropy stability analysis of a SBP-based spatial 
discretization 

Consider a single tensor product element and a spatially discontinuous collocation discretization 
with N = p + 1 solution points in each coordinate direction; the following element-wise matrices 
will be used: 

'D Xl = (T> n 0 In 0 In 0 h) , Dx 2 = (In 0 Dn 0 In 0 h ) , 

Dx 3 = (In 0 In 0 Dn 0 h ) > 

An = (Aat <g> Jw ® /at < 8 > I 5 ) , A X2 = (Ijv 0 Ajy < 8 > -/tv 0 h ) , 

A X3 = (/jv 0 -/at 0 A jv 0 is) , 


T X1 = (Tn ® In ® In ® h) , V X2 = (In ®Vn ® In ® h) , (45) 

Tx 3 = (Cv 0 7at 0 Vn 0 I 5 ) t I’x 1X2 = (Pn 0 Pat is) , 

Pxix 3 = (Pat 0 1a/" 0 0 I 5 ) , Px 2 x 3 = (-Tat 0 TV 0 TV 0 is) , 

P = Pxix 2 x 3 = (TV 0 Pat 0 Pat 0 I 5 ) , 


£> Xl = (£w 0 /w 0 -Tat 0 15 ) 5 £>x 2 = (-Tat 0 Bn ® In ® h ) , 

Sx 3 = (l)v 0 Tv 0 Bn 0 I 5 ) , 

where Pat, TV, Ajv, and Bn are the one-dimensional (ID) SBP operators [40], and Tv is the identity 
matrix of dimension N. I 5 denotes the identity matrix of dimension five. ' The subscripts in (45) 
indicate the coordinate directions to which the operators apply (e.g., T> Xl is the differentiation 
matrix in the x\ direction). Furthermore, we define the norm w'Pq = ||5'||p, where w and S 
are the vector of the entropy variables at the solution points and the mathematical entropy of the 
system, respectively. When applying these operators to the scalar entropy equation in space, a hat 
will be used to differentiate the scalar operator from the full vector operator. For example, 

V = (Vn 0 Pn 0 Vn) • (46) 


Within one tensor product element, the 3D compressible Navier-Stokes equations are discretized 


as 


I + (f</> - rr>) + V£ A, 2 (W> - C) + 


f(i) f(Y) 
r 3 _ r 3 


= v 


-1 

xi 


(B) 

Si 


(In) 

g) 


) + r-(g™ + g<'"> 


+ r 


-1 

x 3 


(B) 

S 3 


(In) 

S 3 


where the vector of conservative variables is ordered as 

q = (? ( x (i)(i)(i)) T . 4 ( s (i)(i)(2)) T , • ■ • , Q (x(N)(N)(N)] T ) = (q(i ) T , Q( 2) T , • • • , 9(at3) T ) , 


(47) 


(48) 


and f and f are the grid fluxes. 7 8 The vector g^ B \ i = 1, 2, 3 enforces the boundary conditions, 

while g 4 - i = 1,2,3 patches interfaces together. The derivatives appearing in the viscous fluxes 

f) J are also computed using the operator T > x . , i = 1, 2, 3 from expression (45). 

7 The 3D compressible Navier-Stokes equations form a system of five non-linear PDEs. 

8 Recall that the vectors with an over-bar are defined at the flux points. 
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As in the continuous case, we apply the entropy analysis to Equation (47) by multiplying 
with w"P from the left. Moreover, we substitute to f ^ , i = 1,2,3 the high-order accurate 
entropy consistent flux constructed according to Equations (33) and (34) with the two-point entropy 
conservative flux presented in Section 4.2.1. The final expression for the time derivative of the 
entropy in the element is then 

-Jl ll'S'll-p + 1 T ^Px2X3^Xl^l + 'PxiX3&X2^2 + 'P X\X2& 3^ 

- w T (V X2X3 B X1 f S' } + V XlX3 B X2 f { p + V X1X2 B X3 4 V) ) + DT (49) 

= w T (v X2X3 ( g[ B) + gS /n) ) + V X1X3 (g ( 2 B) + g^ /n) ) + V X1X2 + g3 /n) )) • 

Note that in (49) the bar over the flux vectors could be safely removed because the contraction of 
(47) against w"P leads only to the boundary fluxes, which satisfy the duality condition (B12) (see 
Figure 1). The quantity DT denotes a positive quadratic term in the first derivative approximation 
of the solution: 

3 3 

DT = ]T (^w) T P[%] (V Xj w) 

i=i j = i 

(V xx w\ T /V[cn] V\ci2\ ‘PfcisA /V Xl w\ 

= I V X2 w I V[c2l] V[C22 \ V[c 23 ] \V X2 W , 

\D X3 W / \P[c 31 ] V[C32 ] V\c 33 }) \D X3 w J 

where [cjj] denotes a block diagonal matrix with blocks corresponding to the viscous coefficients of 
each solution point. The positive semi-definiteness of DT follows from the positivity of the matrices 
Cij used to define [cT,-] (see Appendix B.2 in [16] for the proof). The matrices B Xi , i = 1,2,3 pick 
the interface terms in the respective directions (i.e. , for a high-order accurate scheme on a tensor 
product cell, they pick the solution value at the nodes of the two “opposite” faces). Therefore, 
Equation (49) is the semi-discrete form of Equation (23), which was obtained from the analysis at 
the continuous level. 

5.2 Entropy stability analysis for the solid wall boundary conditions 

We focus now on the construction of an entropy stable penalty term for imposing the solid wall 
boundary conditions for the compressible Navier-Stokes equations. The boundary conditions that 
are developed herein are motivated by the general interface coupling conditions used to cou- 
ple elements in the interior of the domain. The interior conditions combine an entropy stable 
characteristic-based coupling condition for the inviscid terms, with a local discontinuous Galerkin 
(LDG) approach and interior penalty (IP) procedure of the viscous terms. Details of the interior 
coupling terms can be found in Appendix C. These interface terms include a refinement of the IP 
coupling terms relative to those presented elsewhere [41]. 

Without loss of generality, we study a hexahedral element with edge length equal to one and 
we consider only the face plane (0, X 2 ,x 3 ). With these assumptions, Equation (49) reduces to 

J t H-SHp - l Y V X2X3 g {l) ¥ l + w T V X2X3 G {1) f[ V) +DT = w^^pg^. (51) 
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The operators Q^) and G(k) are defined as 


G(k) = { e k ® In ® In) , G(k) = ( e fc ® In ® In ® h ) , (52) 

where 

efc = (0, 0, . . . , 1, 0, . . . , 0,0) T 

is a vector of length N and has a non-zero element corresponding to the solution location, k. 
Therefore, the operators G(k) and G(k) pick out the nodal values of the solution or any flux vector 
at a specific plane according to the ordering introduced in (48). Herein, the face plane (0, X 2 ,xs) 
is characterized by the index k = 1. Thus, Equation (51) represents the contribution to the time 
derivative of the entropy of the boundary points that lie on the face plane ( 0 ,X 2 ,xs). 

In the remainder of this paper, we assume that the node with solution vector q (®(i)(i)(i)) = Q( i) 
(see expression (48)) lies on this face plane. This point will be used to derive the entropy stable 
wall boundary conditions. All numerical states associated to it will be identified with the subscript 

(")(!)’ 

The penalty source term g x is composed of three design-order terms that weakly enforce the 
boundary conditions: 

g[ B > = - (f W (q) - tr (q. g« E >) ) + (f f > - fA B) ) + [M] (w - g™™) . (53) 

In each of the three contributions, the first component (the numerical state) is constructed from 
the numerical solution, while the second component (the boundary state) is constructed from a 
combination of the numerical solution and four independent components of physical boundary 
data. 

The first term enforces the no-penetration Euler wall condition through the inviscid flux of the 
compressible Euler equations. The boundary state is formed by constructing an entropy conserva- 
tive flux based on the numerical state and a manufactured boundary state given by the vector 

9 ^ = 


5 (£) = 


/I 

0 

0 

0 

Vo 


0 0 0 0 \ 

-10 0 0 

1 0 0 

0 1 0 

0 0 1 / 


0 

0 

0 


?5) 


= ( P(l), - (pu l)(l) > {P u 2)(1) , (p«3)(i) 


T 


(54) 


The second term enforces the thermal boundary condition (26), facilitated by manufacturing 
a boundary viscous flux fi’ B . Define the component of the gradient of the entropy variables in 
the numerical state as 0 Xl , 0 X2 and ©^g. Next, specify the value of the heat entropy flux g(t), the 
externally provided bounded function defined as in (26). Finally, define 0 X1 as 

Q X1 = [0^(1), 0 B1 (2), 0^(3), 0^(4), 0 X1 (5)] T , 
where 0^(5) is computed as 

©xi(5) = -g (t) w (1) ( 5) = (55) 

T (i) 
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is constructed as 


With these definitions, the manufactured viscous flux / 


(V,B) 

1 


7(V,B) 
/ 1 


CllQxi T Ci2©x2 “1“ Ci30a; 3 , 


(56) 


where c\j , j = 1,2,3, is constructed from the numerical state. Note that for an adiabatic wall 
g (t) = 0, and from expression (55) we get 0^(5) = 0. 

The third term in (53) enforces the no-slip wall Dirichlet boundary conditions {u\ = 112 = 113 = 
0) through a standard SAT approach. The manufactured boundary state g( NS h Vel is defined in 
terms of entropy variables as 

g (NS),Vei = (^(i), 0, 0, 0, u> (1) (5)) T , (57) 


where io^^I) and 5) are the first and the fifth components of the entropy vector constructed 
from the numerical state. Three boundary conditions are imposed in Equation (57); all velocity 
components are set to zero at the wall. This is immediately clear by recalling that the entropy 
variables for the compressible Navier-Stokes equations are defined as 


w = 



UiUi Ul u 2 U3 

2 ^ ^ rj~i 5 rj~t 5 5 



The matrix [M] is a block diagonal matrix whose five-by-five blocks are defined as 

a< B ) 

M = H cn H, H = diag(l, 1, 1, 1, 0). (58) 

\'xi J (i)(i) 

The matrix cn has the functional form of the usual symmetric positive semi-definite matrix cn 
defined in Appendix A. This matrix has to be constructed using a set of primitive variables that 
is independent of the numerical solution at all times. For instance, for external flows, cn can be 
constructed using the externally provided data at the far-field (e.g., (/?oo 5 |^oo| 5 Koo| 5 1 7?oo I •> Too))- 9 
The coefficient is a positive value used to modify the strength of the SAT penalty term, and 
can be specified by the user. The factor i n the denominator is introduced to achieve 

the correct asymptotic order of accuracy; it allows an increase in the strength of M with increased 
resolution in the normal direction. 

Summarizing Equation (53), the penalty at the face point is the sum of three terms: 

• the difference between inviscid flux and the entropy consistent flux at the node in the normal 
direction; 

• the difference between the internal viscous flux and a boundary viscous flux at the node in 
the normal direction; 

• the difference between the solution (in entropy variables) at the node and the data imposed 
at boundary, multiplied by the matrix M. 

9 In a general framework, the matrix M is built using the five-by-five matrix cn where the index i denotes the 
normal direction to the wall. 
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The penalty term (53) contracted with the entropy variables and simplified, yields the expression 
RHS = - w T V X 2X3 G (1) (i[ 1] (q) - ff (q,g (S) )) 

+ ^ Vx 2 X 3 g {i) (jfJ.fW*)) (59) 

+ W T V X 2 X 3 g { 1 ) [M] (w-g(^ d ). 


The entropy stability of the penalty source term (51) defined by Equation (53) is demonstrated 
in the following theorems. First, the inviscid term is proven to be entropy conservative and then 
entropy stable, if dissipation is added. Next, the second term, which specifies the thermal condition, 
is proven to be bounded by physical data provided by the user. Finally, the third term, which 
specifies the no-slip boundary conditions, is proven to be entropy stable. 

Theorem 5.1. The penalty inviscid flux contribution in Equation (53) is entropy conservative if 
the vector is defined as in (54). 

Proof. The inviscid contribution, to the time derivative of the entropy can be written as (see 
Equations (51) and (59)) 


TN) = (V. 


2 : 22 : 3 ) ( 1 )( 1 ) -^1 W l (^ 2 : 22 : 3 ) ( 1 )( 1 ) / 1 (9(1)) fl (9(1 )i9 


(E) 


(60) 


where W( 1 ) denotes again the entropy variables at the boundary face point and (V X2 x 3 )( 1 )( 1 ) / 0 . 
Substituting the expression for the entropy flux F\ (i.e. , Equation (17) with i = 1) and evaluating 
the entropy consistent flux fff using g (1 ) and yields the desired result 


) — (V x 2 ®s)(l)(l) 

= fPxixfl) ( 1 )( 1 ) 


w T t(7) 
“(i) 


fi’ (9(i)) - Vh - (9i) + 


ST I n niff) 

9(1) 


-ifi + W( 1) f[ r U(i), 9 




= 0. 


(61) 


□ 


Corollary 5.1. The penalty inviscid flux contribution in Equation (53) is entropy stable if the 
vector gi E ) is defined as in (54) and f sr is replaced by the entropy stable flux f ssr defined in (43). 

Remark 5.1. A result similar to Corollary 5.1 is given by Svard and Ozcan [20] in the context of 
high order finite differences approach. 

Using Theorem 5.1 we are left only with the viscous contributions: 

^\\S\\ 2 V + ™ T Vx 2 X 3 G { i)f[ V) + BT < + w t V X 2 X 3 G {1) 

dt \ / ( 62 ) 

+ w T P I «,S (1) [M](w-g< K s>.™). 

Theorem 5.2. The viscous penalty terms in (53) 

are entropy stable for any value of g (t) and any matrix M as defined in (58) . 


21 


Proof. Clearly, the viscous flux term on the left-hand-side (LHS) of (62) is balanced by the same 
term on the RHS. Therefore, expression (62) reduces to 

J t \\S\\ 2 V + DT < - w T V X2X3 g {1 /y B) + W T V X2X3 G(1)[M] (w - g^-™) . (63) 

The contraction — w T 'P X 2 X 3 Cq 1 )f ( j' with defined as in (56) yields the following point-wise 

contribution to the time derivative of the entropy 

— ^( 1 ) (^23;3)(l)(l) /l = (^ 2 : 2 x 3 ) (i)(i) K g(l)i (64) 

when enforcing the heat entropy flow condition (26) through (55). Since g(f) is a known bounded 
function (i.e., L 2 n L°°) expression (64) is also bounded. We highlight that for an adiabatic wall 
g(f) = 0 and consequently the viscous flux penalty in (53) conserves the entropy (as it should) 
because the heat flux is zero. Note that the term (64) mimics exactly the boundary contribution 
that has been obtained from the continuous analysis (see Equation (27)). 

We are then left with the contribution w T Q^[M\ (w — g( NS )’ Vel ^ . At the nodal level, this term 
can be re-written as 


1 - 
2 . 

+ 


w 

1 

2 


T 

( 1 ) 


(fP X 2 X 3 )( 1 )( 1 ) 1 ) 

™ ( 1)-<? (7V5) ’^) T 



^ g (NS),Veiy 

( 1 )( 1 ) M ( W i±) 


(V 1 M n( NS )’ Vel 

\r X2X3J(1)(1) M 9 

g {NS),Vel\ 


(65) 


The penalty contribution given by Equation (65), imposes the no-slip Dirichlet boundary conditions 
on the velocity components, and is bounded if 


• M is negative definite; 


• M is independent of the numerical state. 

If these two conditions are fulfilled, the first and the last term in (65) introduce only dissipation, 
whereas the second one is a bounded term because it is just a function of data, and in general it is 
zero for the no-slip boundary conditions. □ 


For a Reynolds number Re that approaches 00 , we would like to smoothly recover only the 
no-penetration (or wall slip) boundary condition that characterizes the Euler equations (first con- 
tribution in (53)). To achieve that, the five-by-five matrix M needs to be a function of Re and can 
be constructed as in (58), i.e., 

a (B) ~ ,m 

M = — HcuH , H = diag(l, 1, 1, 1, 0), a (B) > 0. 

I 7 Zll(l)(l) 

Recall that cn has the functional form of the usual cn matrix and it is constructed using a state 
that is independent of the numerical solution. 
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6 Numerical results 


The objective of this section is to demonstrate the accuracy and robustness of the new entropy stable 
wall boundary conditions coupled with the family of high-order entropy stable interior operators 
developed in [37]. However, before proceeding with the numerical tests, we demonstrate with 
an example how the application of the standard procedure, which is generally used to imposed 
weakly the solid wall boundary conditions for the linearized Navier-Stokes equations [23], yields 
non-entropy stable boundary treatment. 


6.1 Non-entropy stable no-slip wall boundary condition: Isothermal wall 

In Section 5.2, we have shown that constructing g] as in (53) yields entropy stable wall boundary 
conditions. However, one might attempt to construct g) as the sum of an inviscid penalty flux 
and only a viscous interior penalty term, 


g[ B) = - %) 


f! /} (q) - ff (q, g^)] + 0 (1) [L] (w - g^ S >) 


( 66 ) 


For an isothermal wall, for instance, g^ 5 ) is a vector of data that imposes both the no-slip boundary 
condition on the velocity vector (i.e. , u\ = = 113 = 0) and the wall temperature: 

g {NS) = U (1) (l), 0,0,0,--M . (67) 

\ -Lwall J 

The matrix [L\ in (66) is a block diagonal matrix with blocks of size five-by-five. Comparing the 
two definitions of g x given in Equations (53) and (66), it can be shown that in the latter approach 
no viscous flux penalty terms are introduced. This is a key difference, and as shown in Appendix 
D, yields non-entropy stable solid wall boundary conditions that are highly unstable when used in 
combination with fine grids and/or high-order accurate polynomial representations of the solution. 


6.2 Computation of a square cylinder in subsonic flow 

The flow past a square cylinder represents a benchmark test case for external flow past bluff 
bodies. This flow has been the subject of intense experimental and numerical research in the past. 
In fact, this bluff body is a simple but a central shape for many engineering applications, including 
aeroacoustics and air pollutant transport and dispersion in urban environments. 

The flow is described in a Cartesian coordinate system (aq, X 2 , X3), in which the aq-axis is 
aligned with the inlet flow direction, the aq-axis is parallel with the cylinder axis and the aq-axis 
is perpendicular to both directions (see Figure 2). A fixed two-dimensional square cylinder with 
a side d is exposed to a uniform freestream velocity vector with modulus rioo- The length of the 
square cylinder in the aq-direction is 10 d. 

The following boundary conditions are used. A uniform flow is prescribed at the inlet which 
is located 10 d units upstream of the cylinder. At the outlet, located 20 d unit downstream of the 
cylinder, farfield boundary conditions are used. A no-penetration (Euler) boundary condition is 
prescribed at the upper and lower boundaries. No-slip and adiabatic conditions are enforced at 
the body surface. A periodic boundary condition is used in the spanwise direction aq. In the 
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^-direction, the solid blockage of the confined flow (i.e. , the vertical distance between the upper 
and the lower inviscid walls) is set to 18 d. 

The flow has a freestream Mach number of M 0 0 = 0.1, and a Reynolds number of Re ^ = 2 x 10 2 . 
The Reynolds number is based on the modulus of the freestream velocity vector Uoo and the height 
of the cylinder d. At this Reynolds number, the regime is laminar and it usually persists up to 
a Reynolds number of about 4 x 10 2 . Moreover, the vortex shedding is characterized by one very 
well-defined frequency [42]. A very small time step is used to integrate the system of ordinary 
differential equations (ODEs) so that the temporal error is negligible compared to that of the 
spatial discretization. 

6.2.1 Accuracy of the no-slip wall boundary conditions 

The proposed entropy stable no-slip wall boundary conditions do not force the numerical solution 
to exactly fulfill the boundary conditions. Instead the effect can be described as a rubber-band 
pulling the solution towards the boundary conditions. The computed boundary value (or numerical 
state) typically deviates slightly from the prescribed value but the deviation is reduced as the grid 
is refined. Therefore, the error at the boundary can serve as a rough measure of the error of the 
entire solution. 

We compute the maximum norm L°° of the error of the three velocity components u\ , U2 and 
U3 on the complete surface of the cylinder at t = 1, for three different grids. The meshes are fully 
unstructured, although a structured subdivision is used around the square cylinder and the near 
wake region to perform a grid convergence study (see Figure 2). Grid 3 has 20 points on each side 
of the square, 20 points in the “radial” direction in structured portion near the body, 40 points in 
the near wake region in the freestream direction, and 8 points in the spanwise direction. Grid 2 and 
Grid 1 are obtained by taking every other and every fourth grid point of Grid 3 in the structured 
region. The simulations are performed using different orders of the polynomial ( p = 1,2, 3, 4). The 
results are shown in Tables 1, 2 and 3. 



Figure 2. Example of structured/unstructured grids use for the flow past a 3D square cylinder at 
Re oo = 2 x 10 2 and M ^ = 0.1. 

We highlight a few observations. First, in all cases an increase in theoretical order of accuracy 
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results in a error reduction on all grids. Secondly, although the convergence rates in model problems 
are shown on much finer meshes, the computed order of accuracy is very close to the formal value 
between the medium and fine grids, even for these more realistic meshes. 


Table 1 L°° error norm of the velocity component u\ at the wall and convergence rates; t = 1; 3D 
unsteady laminar flow past a square cylinder at Re oo = 2 x 10 2 and = 0.1. 



p = 1 

rate 

P = 2 

rate 

p = 3 

rate 

p = 4 

rate 

Grid 1 

4.73e-2 

- 

2.15e-2 

- 

9.61e-3 

- 

4.54e-3 

- 

Grid 2 

1.47e-2 

1.69 

2.88e-3 

2.90 

5.83e-4 

4.04 

1.39e-4 

5.02 

Grid 3 

3.55e-3 

2.04 

3.66e-4 

2.98 

3.40e-5 

4.10 

4.52e-6 

4.94 


Table 2 L°° error norm of the velocity component U 2 at the wall and convergence rates; t = 1; 3D 
unsteady laminar flow past a square cylinder at Re 0 c = 2 x 10 2 and = 0.1. 



p = 1 

rate 

p = 2 

rate 

p = 3 

rate 

p = 4 

rate 

Grid 1 

7.20e-2 

- 

2.71e-2 

- 

1.43e-2 

- 

5.14e-3 

- 

Grid 2 

1.79e-2 

2.01 

3.34e-3 

3.02 

1.10e-3 

3.70 

1.96e-4 

4.71 

Grid 3 

4.65e-3 

1.94 

4.20e-4 

2.99 

7.21e-5 

3.93 

6.48e-6 

4.92 


Table 3 L°° error norm of the velocity component 113 at the wall and convergence rates; t = 1; 3D 
unsteady laminar flow past a square cylinder at Re ^ = 2 x 10 2 and = 0.1. 



p = 1 

rate 

P = 2 

rate 

p = 3 

rate 

p = 4 

rate 

Grid 1 

2.75e-4 

- 

1.34e-4 

- 

1. Ole-4 

- 

8.62e-5 

- 

Grid 2 

5.98e-5 

2.20 

1.71e-5 

2.97 

7.92e-6 

3.67 

3.14e-6 

4.78 

Grid 3 

1.38e-5 

2.12 

2.03e-6 

3.04 

5.30e-7 

3.90 

8.70e-8 

5.17 


6.2.2 Vortex shedding 

In this section we investigate the vortex shedding and the time variation of the lift and drag 
coefficients. We compare our results against the data reported by Sohankar et al. [43]. We compute 
the following quantities: The Strouhal number, f d/uoo, where / is the frequency of the vortex 
shedding; the tinre-averaged drag coefficient, cd, and the root-nrean-square (RMS) of the spanwise- 
averaged lift coefficient, c^ MS . We use the same grids presented in the previous section, and 
different orders of the polynomial ( p = 1,2, 3, 4). The results are illustrated in Tables 4, 5, 6. From 
these tables, it can be seen that in all cases the accuracy of the results improve by increasing the 
order of accuracy of the scheme. We also note that, on Grid 3, which is very coarse compared 
to the typical grids used with second-order FV and FD schemes, fourth- (p = 3) and fifth-order 
(p = 4) accurate entropy stable schemes perform very well. In fact, the aerodynamic coefficients 
computed with these two discretizations are in very good agreement with the results reported in 
literature [43]. 
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Table 4 Strouhal number, mean drag coefficient and spanwise-averaged RMS of the lift coefficient 
for the 3D unsteady laminar flow past a square cylinder at Re oo = 2 x 10 2 and Moo = 0.1; Grid 1. 


Solution 

St 

( CD ) 

n RMS 

C L 

SSDC p = 1 

0.098 

1.01 

0.02 

SSDC p = 2 

0.109 

1.08 

0.06 

SSDC p = 3 

0.142 

1.19 

0.11 

SSDC p = 4 

0.151 

1.28 

0.15 

Sohankar et al. [43] 

0.160 

1.41 

0.22 


Table 5 Strouhal number, mean drag coefficient and spanwise-averaged RMS of the lift coefficient 
for the 3D unsteady laminar flow past a square cylinder at Re oo = 2 x 10 2 and Moo = 0.1; Grid 2. 


Solution 

St 

( C D ) 

„RMS 

C L 

SSDC p = 1 

0.128 

1.16 

0.07 

SSDC p = 2 

0.139 

1.28 

0.13 

SSDC p = 3 

0.153 

1.36 

0.20 

SSDC p = 4 

0.159 

1.40 

0.23 

Sohankar et al. [43] 

0.160 

1.41 

0.22 


Table 6 Strouhal number, mean drag coefficient and spanwise-averaged RMS of the lift coefficient 
for the 3D unsteady laminar flow past a square cylinder at Re oo = 2 x 10 2 and Moo = 0.1; Grid 3. 


Solution 

St 

( c d) 

n RMS 

C L 

SSDC p = 1 

0.134 

1.29 

0.12 

SSDC p = 2 

0.154 

1.37 

0.19 

SSDC p = 3 

0.159 

1.40 

0.22 

SSDC p = 4 

0.159 

1.42 

0.23 

Sohankar et al. [43] 

0.160 

1.41 

0.22 
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6.3 Computation of a square cylinder in supersonic free stream 

The development of a high-order accurate entropy stable discretization aims to provide the next 
generation of robust high fidelity numerical solvers for complex fluid flow simulations, for which 
standard suboptimal algorithms suffer greatly or fail completely. By computing the flow past a 3D 
square cylinder at Re oo = 10 4 and M m = 1.5, we provide numerical evidence of such robustness for 
the complete entropy stable high order spatial discretization. This supersonic flow is characterized 
by a very large range of length scales, strong shocks and expansion regions that interact with each 
other, leading to complex flow patterns. During the past three decades, this fluid flow problem has 
been thoroughly investigated by several researchers for aerodynamic applications (see for instance 
[44-46]). 

The domain of interest spans one square cylinder edge in the X 3 direction, and at the two 
planes perpendicular to this coordinate direction, periodic boundary conditions are used. The 
flow is computed using an unstructured grids with 43, 936 hexahedrons. A fourth-order accurate 
(p = 3) entropy stable discretization without any stabilization technique is used. The body surface 
is considered adiabatic. The solution is initialized using a uniform flow at M ^ = 1.5 with zero 
angle of attack. 

At the beginning of the simulation a strong shock is formed in front of the blunt body. Subse- 
quently, the discontinuity moves upstream until it reaches a stationary position that is about 2.15 
square cylinder edges far from the frontal surface of the body. During this phase, additional weaker 
shocks, which originate from the four sharp corners of the blunt body, interact with the subsonic 
regions formed near the walls. This complicated flow pattern, yields the formation of shock-lets in 
the wake of the square cylinder. Figure 3 shows the high order grid near the blunt body and the 
Mach number and density contours at t = 1.5. It can be seen that relatively small oscillations are 
generated in front of the shock. This numerical feature is absolutely natural and expected because 
the solution has been computed with a fourth-order accurate scheme without artificial dissipation 
or filtering technique. Nevertheless, the simulation remains stable at all time, and the oscillations 
are always confined in small regions close to the discontinuities. 

In Figures 4 a global view of the high order grid, the Mach number, density, temperature and 
entropy contours at t = 100 are shown. At t = 100, the shock has already reached a stationary 
position, and the flow past the square cylinder is completely unsteady, characterized by subsonic 
and supersonic regions. The formation of shock-lets in the near wake region are clearly visible. 


7 Conclusions 

Herein, we have shown that a no-slip boundary condition together with a boundary condition on the 
heat entropy flow, (1/T dT/dn) wall , imply stability for the continuous compressible Navier-Stokes 
equations. The boundary condition on the heat entropy flow is in complete agreement with the 
thermodynamic (entropy) analysis of a generic system. An entropy stable numerical procedure is 
presented for weakly enforcing these solid wall boundary conditions via a penalty approach. The 
resulting semi-discrete operator mimics exactly the behavior at the continuous level. The proposed 
non-linear boundary treatment provides one of the most important missing information for com- 
pleting the non-linear stability in the L 2 norm of the continuous and semi-discretized compressible 
Navier-Stokes equations. Although discontinuous spectral collocation operators are used in this 
work, the new boundary conditions are compatible with any diagonal norm summation-by-parts 


27 



(a) High order grid in the near-body and near-wake. 



(b) Mach number; AM = 0.0146. 



(c) Density; A p = 0.0114. 


Figure 3. Unsteady flow past a 3D square cylinder at Re oo = 10 4 and = 1.5, computed with 
the fourth-order ( p = 3) accurate entropy stable spatial discretization; t = 1.5. 
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(a) High order grid. 



(b) Mach number; AM = 0.0095. 



(c) Density; A p = 0.0090. 
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(d) Temperature; AT = 0.0077. 



(e) Entropy; As = 0.0044. 


Figure 4. Unsteady flow past a 3D square cylinder at Re oo = 10 4 and = 1.5, computed with 
the fourth-order (p = 3) accurate entropy stable spatial discretization; t = 100. 
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spatial operator, including finite element, finite volume, finite difference, discontinuous Galerkin, 
and flux reconstruction schemes. 

Numerical computations around a three-dimensional square cylinder in the subsonic regime are 
performed to highlight the accuracy and robustness of the proposed numerical procedure. Mea- 
surement of forces on the cylinder showed very good agreement with the results available from the 
literature. Furthermore, we have shown that wall penetration velocity (admissible in a penalty 
technique) approaches zero to design-order. 

The robustness of the complete semi-discrete operator (i.e. , the entropy stable interior operator 
coupled with new boundary conditions) has been demonstrated for the supersonic flow past a three- 
dimensional square cylinder at Re oo = 10 4 and = 1.5. This test has been successfully computed 
with a fourth-order accurate method without the need to introduce artificial dissipation, limiting 
techniques or filtering, for stabilizing the computations, a feat unattainable with several alternative 
approaches to wall boundary conditions and high order interior operators based on linear analysis. 

This work clearly indicates that, although incremental improvements to existing algorithms 
will continue to improve overall capabilities, the development of novel robust numerical techniques 
such as entropy preserving or entropy stable schemes and their extension to complex problems of 
industrial relevance offers the possibility of radical advances in this area in terms of robustness, 
fidelity and efficiency. 
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Appendix A 


Coefficient matrices of the viscous flux 


The viscous coefficient matrices A used to define the viscous fluxes in Cartesian coordinates in 

L J 

(7) are defined as 
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The symmetrized coefficient matrices used in (13) to define the viscous fluxes as a function of 
the gradient of the entropy variables are found using [16] 
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Appendix B 


Summation-by-parts operators 


With most methods, the discretization of the domain P is obtained by a space-filling non- 
overlapping subdivision into small volumes, called cells or elements. Such a subdivision is called 
a grid or mesh. In this section, without loss of generality, we introduce the summation-by-parts 
operators for a generic ID spatial discretization and specialized them to high-order accurate dis- 
continuous collocation methods. As shown in the next sections, tensor product algebra allows the 
results to extend directly to three-dimensions. 

Define on the interval — 1 < x < 1, the vectors of N discrete solution points 

X = ■ ■ ■ ,x N -i,x N ) T , -1 < Xi < X 2 , ■ ■ ■ , < X N -1,< X N < 1 • (Bl) 

Because the approximate solution is constructed at these points, they are denoted as solution points. 

B.l First derivative 

First derivative operators that satisfy the SBP convention, discretely mimic the integration-by-parts 
condition 

xr xr 

(B2) 

XL X L 

where 4> = 4> (x) is a generic test function of the independent variable x. This mimetic property is 
achieved by constructing the first derivative approximation of 4> at IV solution points, 

V4>, = ),(j)(x 2 ),...,(j){xN)) T , 

with an operator in the form 

V = V~ 1 Q 1 v = v T , c T vc>0, c/o, 

Q T = B-Q, B = diag (—1,0, ... ,0, 1) . 

While it is not true in general that V is diagonal, herein the focus is exclusively on diagonal norm 
SBP operators, based on fixed element-based polynomials of order p. The matrix V incorporates 
the local grid spacing into the derivative definition and it may be thought of as a mass matrix in 
the context of Galerkin FEM. The nearly skew-symmetric matrix, Q, is an undivided differencing 
operator where all rows sum to zero and the first and last column sum to —1 and 1, respectively. 

Integration in the approximation space is conducted using an inner product with the appropriate 
integration weights contained in the norm V, 

XR 

J fi^dx & <f) T VVci, q= (g(xi),g(x 2 ),...,g(xAr)) T . (B4) 

XL 

Using the definition in (B3), the SBP property can be demonstrated: 

< t^VV-'Qq . = cf) T (b - Q T ) q = (j> N q N - Ml ~ 4> T V T V q. (B5) 

The specific first derivative SBP operator used in this work is shown in Appendix B.5.1. 
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B.2 Second derivative 


The viscous approximations also satisfy the SBP condition. Integration by parts yields 




XR 


XL 



dq 

dx 


dx, 


XL 


(B6) 


where id = id(x). The second derivative variable coefficient operator resulting from two applications 
of the first derivative may be manipulated for diagonal norm, V, into the expression 

V 2 {d) = V~ v [~V T Vi JV + BfiV} , V T VfiV = (v J VtiVy , 

•& = diag ($(x)) , (B7) 

C T (v T V'dv^ c > o, c T &C > o, vc. 

The P-norm inner product yields the expression 

0 T 7TP _1 (-V'T&V + BtfV^J q = 4> T BtiV q - 0 T (v T VtiV} q, (B8) 

which is usually the form used to prove stability of the viscous terms. It is clear that the continuous 
interface terms are mimicked. Likewise, based on the definition (B7), the second term on the RHS 
of Equation (B6) can be approximate as 


/ 


dx^^dq 
dx dx 


dx « 


XL 


0 T (t^WV ) q. 


B.3 Complementary grids 

The entropy analysis that will follow is performed on a staggered set of solution and flux points. 
The flux points a set of (N + 1) intermediate points prescribing bounding control volumes about 
each solution point (see Figure 1). These points are similar in nature to the control volume edges 
employed in the finite volume method. The distribution of the flux points depends on the dis- 
cretization operator. The spacing between the flux points is implicitly contained in the norm V. 
In fact, the diagonal elements of V are equal to the spacing between flux points, 

X = (x 0 ,Xi,...X7v) T , Xo = Xi, XJV = XN, 

(By] 

Xi Xi— 1 ^(i)(i) ; ^ f? ■ ■ • , 


In operator notation, this is equivalent to 


Ax = VI, 


(BIO) 


where 

1 = ( 1 , 1 ,..., 1 ) T , 
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is a vector with N elements and 


/ —1 1 0 

0 -1 1 


0 0 0 

\ o o o 


0 0 0 \ 

0 0 0 

0 0 

-1 1 0 

0 - 11 / 


(Bll) 


is used to calculate the undivided difference of the two adjacent fluxes. All the quantities located 
at flux points are denoted with an over-bar. Note that in (B9), the first and last flux points are 
coincident with the first and last solution points, which enables the endpoint fluxes to be consistent, 
i.e., 

/o = f(qi), In = /(<?a0- (B12) 

This duality is needed to define unique operators and is important in proving entropy stability [37]. 


B.4 Telescopic flux form 

All SBP derivative operators, T>, can be manipulated into the telescopic flux form, 

=v- 1 Q{ + T (p+1) = V-'Ai + T (P+1) , (B13) 

where T( v +i) is the truncation error of the approximation. This form admits a generalized SBP 
property: All SBP operators defined in Equation (B3) can be manipulated to transfer the action 
of the discrete derivative onto a test function with an equivalent order of approximation [47]. 

Likewise, the variable coefficient viscous operators presented in Section B.2 may be expressed 
in the form 

— ~ p-i (-V T VVV + BtiV) q = p-' Af (,)) , (B14) 

ax \ ax J V / 

and satisfy a telescoping conservation property, which is identical to that of the inviscid terms in 
(B13). 


B.5 Spectral discretization operators 

Spectral collocation methods are commonly implemented on computational grids based on the nodes 
of Gauss-quadrature formulas (i.e., Gauss, Gauss Radau, or Gauss Lobatto (GL)). These smooth 
but nonuniform grids are highly clustered at the boundaries of the domain, in stark contrast to the 
uniform grids used in conventional finite difference methods. 

The numerical methods developed herein are all collocated at the Legendre-Gauss-Lobatto 
(LGL) points, and include both end points of the interval. This distribution which includes the end 
points, allows the operators to be written in terms of flux differences, analogous to a finite volume 
method and consistent with Equations (B13) and (B14). The complete discretization operator for 
the p = 4 element is illustrated in Figure 1. 
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B.5.1 Lagrange polynomials 

Define the Lagrange polynomials on the one-dimensional N discrete points x as 


' L(x) = (x- x 2 )(x - x 3 )...(x - x n - 2 )(x - X N -l) 
(1 — x)L(x) 


L\(x) = 
Ln(x ) = 
Lj{x) = 


2Z(-1) 

(1 + x)L(x) 

2L(+1) 

(1 — x 2 )L(x) 


2<j<N-l. 


(x — Xj)L'(xj) 

With a slight abuse of notation, define the vector of Lagrange polynomials as 


L(x) = (Li(x),L 2 (x), • • • ,Lat_i(x),L7v(x)) t . 


(B15) 


(B16) 


A Legendre collocation operator may be constructed by approximating integrals by the LGL quadra- 
ture formula (see [37] for a detailed derivation). Let rj = (771 , rj 2 , ■ ■ ■ r]N-ii w) T be the nodes of 
the LGL quadrature formula (i.e. , the zeroes of the polynomial P! n _ l {x){ 1 — x 2 ) [48]), and let 
u ;j, 1 < l < N be the quadrature weights. 

The mass and stiffness matrices V and Q are defined as 

V = y^L(? 7 ;;x) [L(rji;x)] T u;i, (B17a) 

1 

Q = ^2Mvr,x)[L'(rn]x)] T ui, (B17b) 

1 

where L '(j]i\x) denotes the derivative of the Lagrange polynomials with respect to x. The matrix V 
is SPD for any x. Given these definitions, the differentiation matrix V is then computed as usual: 

V = V~ l Q. 


The discrete operators are provided for polynomial orders one to four in Table Bl. Only the 
upper triangular portion of Q is provided. The full Q matrix may be reconstructed from the 
skew-symmetry property Q + Q 1 = B. 
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Table B1 Differentiation operators for polynomials of degree p one through four. 
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Appendix C 


Entropy stable discontinuous interfaces 


In this Appendix we focus on the development of entropy stable interface penalties to patch 
together two elements. The interior conditions combine an entropy stable characteristic-based 
coupling condition for the inviscid terms, with a LDG approach and IP procedure for the viscous 
terms. Consider two hexahedral elements whose faces are all internal faces. Thus, for these two 
elements there is no boundary contribution to the entropy estimate. Then, Equation (47) yields 


fW p i p— 1 A f(-0 T) f (V) 

-fa + '' xi,l A zi,dl,Z - D xi Pu + X 2 ,l ~ D X2,l t 2,l 

+ A",!i a «.Aw 

— ^ 1 x 2 ,l&2,l ^ ' x 3 ,l&3,l > 

^ r i p-i a _ p f ^ ) i p— i a f ^ _ p 

' ' xi,r L -^xi,r 1 l,r ,r 1 l,r w / X2i? . ‘-±x 2 ,r 1 2,r 1 2,r 


(Cla) 


i—l A <?C0 -n f?00 
1 3,r 


T ^x 3 ,r ^3 } r 

= p-l „(*"). p-1 G™) , p-1 e ( /n ) 

' xi,r &l,r 1 ' X2,r o 2 ,r 1 ' xs,r 03^ j 


(Clb) 


where the subscripts l and r denote the left and right elements. The coupling between the two 
cells is achieved through the penalty interface terms gj / / "' 1 and g^”\ which are constructed as a 
combination of an LDG-type approach and an IP technique. Therefore, Equations (Cl) can be 
re-written as 
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(C2c) 


(C2d) 


where ©p and @j r are the vectors of the gradient of the entropy variables on the left and right 
elements in the i direction, and are the penalty interface terms on the variable q. and 

the gradient of the variable ©., respectively. The matrices [cip] are the block diagonal coefficient 
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matrices that have been introduced in Equation (50). Note that in (C2) the Einstein summation 
convention notation is used. 

Multiplying the two discrete equations in the left element by w J'Pi and ([%,;] respec- 

tively, and the two discrete equations in the right element by w ( T 'P r and ([%,r]©j,r) T 'Pr, respec- 
tively, adding to each equations its transpose and using the SBP property Qn = \Bn + Skew (Qn), 
the final expression for the time derivative of the entropy S in each each element is 


s' + 2 

1 (^Px2X3,l T^X±XS,l 13x2,1^2,1 T^XlX2,l ^X3,l^3,l^ 

= W (VX 2 X 3 ,1 B £i,z[ c l j,l] ®j,l *b T*xi x 3 ,l ^X2,l[ c 2j,l] ®j,l d" 'P X 1 0C2,1 B X3 ,l [ c 3 j,i] ©7,2 ) 
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(C3a) 


dt ll5rll ^ + 2 


\] [ °ij,r \ ®j,r 
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(C3b) 

To simplify the notation (but without loss of generality), assume that each face of the elements 
is parallel to one of the axes of a Cartesian reference system and the interface between the two 
cubes lies at x\ = 0. We also assume that all the points that lie on the other faces of the two cubes 
are treated in an entropy stable fashion. Then, for our analysis we can just focus on the interface at 
x\ = 0. Therefore, introducing again the operators and G(k)> where k = N for the left element 
and k = 1 for the right one, one arrives at 
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The interface penalty terms are constructed as follow: 
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(C5b) 


(C5c) 


(C5d) 


The superscripts (— ) and (+) denote the collocated values on the left and right side of the interface, 
respectively. The LDG penalty terms involve the coefficients ^(l±a) and act only in the normal 
direction to the face. The IP-FEM terms involve the parameter block diagonal matrix, [y], whose 
five-by-five block matrices are left unspecified for the moment. Substituting expressions (C5) in 
(C3) and summing all the contributions of the two elements in a single equation results in 

U7 ll'S'/ll'P; + ll'Sr||-p r +2 J [ Cij ,i\ ®j.l + J [Cij : r] &j,r ^ + T^, (C6) 

az i L Vl Vr_ 

where Y ( + and Y ( ' ' are the inviscid and the viscous interface terms. At each node, these terms 
Y+ } = [w^ 1 — f ssr Qi + ^ — ^ 


are 


= (w^ — w ^ -* 


^ A — vj( ^ , 

T (1 } = (w^ - Y (w^ - w- +) ) , 


(C7) 

(C8) 


where A and Y are five- by- five matrices. Clearly, the interface contributions are dissipative if both 
A and Y are negative definite. The matrix A is constructed as in (43). As for the no-slip boundary 
conditions case, we would like the matrix Y to go to zero for Re —> oo. To achieve that, the matrix 
Y needs to be a function of the inverse of the Reynolds number. In this work, we construct this 
matrix as follows 

y = -a (/n) ,^L, + . 11 , (C9) 

2 (7 7 a: 1 )( 1 )( 1 ) 

where c^ 1 ' ) and are the matrices at the left and the right side of the interface in the normal 
direction. The coefficient a^ In ^ is a positive value that can be used to modify the strength of the 
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IP penalty term. The factor (7^)/-^-^ in the denominator has been introduced to get the correct 
dimension and it increases the strength of N with increased resolution in the normal direction. 

Remark C.l. The parameter values a = 0 and a = ±1 yield a symmetric LDG and a “ flip- 
flop ” narrow stencil LDG penalty, respectively. An LDG value of a = 0 produces a global discrete 
operator that has a neutrally damped spurious eigenmode. The IP dissipation effectively damps this 
mode. 
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Appendix D 


Counter example of non-linear wall boundary conditions 


Carrying out the entropy stability analysis by using expression (66), it can be easily shown that 
the inviscid penalty in (66) is entropy conservative (see the proof of Theorem 5.1). Therefore, the 
remaining relation to analyze is 

j t \\Sf v + w T 'P X2X3 G(i)fi' > + DT < + w T V X2 x 3 G (i)[L] — g( NS ^ . (Dl) 

To obtain a quadratic form in boundary terms we need to borrow from DT (see Equation (50)): 

3 3 

DT = E E ^> W ) T w ) 

i = 1 j = 1 

/V xl w\ T /V[cn] V[c 12 ] V[c 13 }\ fV Xl w\ 

= P X2 W V[c 2 l] 'P[c 22 j V[C 2 3 ] V X2 W 

\D X3 wJ \7^[c3l] V[c 3 2] ^[ 033 ]/ \V X3 w/ 

/®ii w \ /T’o;23;3[cil] 'P X2 x 3 [ci 2 ] Px 2 X3 [ c 13]\ /Tij w\ 

= £> X2 W (T’jji )(!)(!) Pxzxs [c 2 l] T , X 2X 3 [C23] T> X2 w +DT 

\Dx 3 w / \Px2X3 [ c 3l] Px 2 x 3 [ c 32] ^2X3 M/ \^x 3 w / 

/P xl w\ T /[cn] [C 12 ] [cis]\ /T> Xl w\ _____ 

= T» X2 w [c 2 i] [c 22 ] [c 23 ] T» X2 w + DT, 

w/ \[c3i] [c 32 ] [C 33 ]/ \T> X3 w / 

where 

7^ = diag (T > X2 X 3 iPx 2 x 3 iPx 2 x 3 ) j 

and the scalar ('P Xl )( 1 ^ 1 ) > 0. Therefore, Equation (Dl) may be written as 


> IH- +DT 


< 


1 ~T-fy 

2 w p 


+ 5" ,T r«i s 6(i)[r] w - w 1 P* 2 x,§(i)[i] g 


W -pn] -P 12 ] 

[ c n] (T’ Xl )( 1 )( 1 ) [cn] ^2 ('P Xl )^ 1 ^ 1 ^ [ci 2 ] 
-[C 12 ] -2(P X1 )( 1 )( 1 ) [c i2 ] -2(7 ? X1 ) (1)(1 ) [c 22 ] 
[ c i 3 ] ~ 2 (T’x!)^)^) [C 13 ] ~2 ('P X1 )^ 1 ^ 1 ^ [c 23 ] 


T -1 


r (JVS) 


-[C 13 ] \ 

— 2 (T’xi)^)^) [ci 3 ] 

— 2 (Tln)^^!) [c 23 ] 

— 2 ("Px! )(!)(!) [C33] y 


w 


(D3) 


where 


and 


T >/ = diag (V X2X3 ! T* x 2 X3 ) 'P X 2 X3 > P X 2 X3) 1 


w T = 


w 


„T 

(0,x 2 ,x 3 ) 


,0,(P 


W)m , (V X2 W) J 0 ^ X3) , (V X3 W) J i x 2i x 3 ) , 0 


w qo,x2,x 3 ) 


(D4) 
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The bold zeros in (D4) indicates that all the numerical states (entropy variables and gradients of 
the entropy variables) of the nodes which do not lie on the plane (0, X 2 ,xs) are set to zero. 

To bound the time derivative of the entropy we must ensure that each term is bounded. The 
first contribution on the RHS is a quadratic term in w and dissipative if the large matrix in (D3) is 
symmetric negative semi-definite. However, to ensure that we only need to construct the following 
20 x 20 matrix, 



( L 

-Cll 

Cl 2 

-C 13 ^ 

1 

-Cll 

—2 (V Xl ) cn 

2 (Vx i)(i)(!) Cl 2 

— 2 (V X \ ) (1)(1) C 1 3 

2 

Cl 2 

—2 (V X 1 ) (i)(i) c i 2 

2 (P^ )(i)(i) C22 

2 C^xi )(i)(i) C 23 


y-Cl 3 

— 2 ('P Xl )(1)(1) c 13 

— 2 {V Xl )(!)(!) c 23 

— 2 (V x i )(i)(i) C33 J 


(D5) 


so that it is symmetric negative semi-definite. The matrices Cij,i,j = 1,2, 3, in (D5) are constructed 
using the primitive variable at the usual boundary node. The rows and columns of the matrix T 
corresponding to the density components are all zero because the first component of the viscous 
fluxes is zero. Therefore, such rows and columns do not affect the negativity of (??). The matrix 
r can be expressed in block form as 

r =(fiT d )' (D6) 

The condition on the five-by-five matrix L that ensures the negative-definiteness of (D6) can ob- 
tained by requiring that the Schur complement of T is negative, 


Schur = A- BD~ l B T = 


cn + 2L (V Xl ) m) ^ 

^ (v xi ) im) 


(D7) 


The inequality in (D7) is a sufficient condition because the block matrix D is already well behaved 
(i.e., it is already symmetric and positive semi-definite). Thus, the Schur complement (D7) is 
smaller than or equal to zero if 


L < 


Cll 

2 ('Pxi )(!)(!) 


(D8) 


Expression (D8) is very similar to relation (66), but it yields a non-entropy stable boundary treat- 
ment, as shown below. 

The last two terms on the RHS of expression (D3) are the remaining contributions to bound. 
Such terms can be re-written in a quadratic form as 


7,™ T/ P X2X3 G(i)[L\ w - w T V X2X3 G {l) [L\ g (A,S) = + \ ( w - g (7V5) ) V X2X3 G (1) [L] (w - g (A " S) ) 

-^(g (iYS) ) T ^3^)Wg (Ar5) - 

(D9) 

From inequality (D8), we know that the matrix [L] is negative semi-definite. Therefore, the first 
term, which is quadratic in (w — g (Ars ^), is dissipative. However, the second contribution is a 
positive term and cannot be bounded because it is not only a function of the imposed boundary 
data In fact, the element in the fifth row and fifth column of the matrix L is non-zero and 

it is a function of the numerical solution through relation (D8) (i.e., through the matrix cR which 
is built from the numerical state at the boundary node). 
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